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ABSTRACT 


A  shell  formulation  was  developed  from  a  three-dimensional  solid.  The  shell  ele¬ 
ment  is  an  isoparametric  element,  and  has  four  corner  nodes  at  which  there  are  three 
displacements  and  three  rotations,  independently.  Therefore,  the  element  formulation 
includes  the  transverse  shear  deformation  and  the  transverse  normal  deformation.  In 
addition,  the  formulation  consists  of  separate  components  of  the  mean  stress  and  devi- 
atoric  stresses  because  the  Gurson  constitutive  model  for  void  growth  is  based  on  the 
mean  stress  and  the  dilatation.  As  a  result,  the  Gurson  void  model  can  be  implemented 
in  the  shell  formulation  at  the  next  stage.  The  shell  element  uses  the  reduced  integration 
along  the  inplane  axes  and  full  integration  along  the  transverse  direction.  If  more  accu¬ 
racy  is  required  along  the  thickness  of  the  shell,  a  large  number  of  integration  points  can 
be  selected  in  the  direction.  Verification  of  the  shell  element  was  performed  for  a  plate 
problem  and  a  shell  problem  whose  analytical  solutions  were  available.  The  next  phase 
of  the  work  is  to  implement  the  Gurson  model  into  the  shell  element.  Once  those  are 
successfully  completed,  the  module  will  be  incorporated  into  the  DYSMAS  program. 
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1.  INTRODUCTION 


Computer  simulation  of  ship  shock  trial  is  one  of  the  goals  in  the  U.S.  Navy  even  if 
there  are  some  Hxiutations  in  the  simulation  because  some  phenomena  cannot  be  repre¬ 
sented  by  mathematical  models.  However,  replacing  physical  testing  by  computer  simula¬ 
tion  as  many  as  possible  will  reduce  the  costs  in  the  design,  production,  and  quahfication 
processes. 

As  a  result,  a  computer  simulation  program  has  been  developed  by  the  Naval  Surface 
^Varfare  Center  (NSW^C)  in  collaboration  with  Germany.  The  program  is  expected  to 
model  both  linear  and  nonhnear  behaviors  of  ship  structures  under  a  shock  loading.  In 
other  words,  failure  process  needs  to  be  modeled  in  a  accurate  and  reliable  way.  One 
of  the  efforts  to  model  such  a  nonhnear  failure  process  is  to  implement  the  damage 
constitutive  equations  in  the  program  such  as  Gurson’s  void  model. 

Because  the  major  ship  structure  is  a  shell,  shell  elements  are  used  for  modeling.  Fur¬ 
ther,  as  stated  above,  the  shell  element  must  include  the  dmamage  constitutive  equation 
like  Gurson’s  void  model.  Therefore,  the  overall  objective  of  this  research  is  to  develop 
a  shell  element  which  can  include  the  damage  constitutive  equation. 

The  research  consists  of  a  couple  of  different  phases.  The  first  phase  is  to  develop 
a  shell  element  which  can  include  Gurson’s  void  model  at  the  next  phase.  One  of  the 
requirements  for  the  shell  formulation  is  that  the  formulation  should  be  able  to  represent 
the  pressure  (mean  stress)  and  deviatoric  stresses  separately  because  Gurson’s  void  model 
was  based  on  the  pressure  component. 

The  second  phase  is  to  implement  the  void  model  into  the  shell  formulation  developed 
in  the  first  phase.  Finally,  the  last  phase  is  to  incorporate  the  whole  module  into  the 
DYSMAS  program.  For  each  phase  of  work,  extensive  verifications  of  the  developed  shell 
element  must  be  conducted.  This  report  corvers  the  work  for  the  first  phase. 
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2.  FINITE  ELEMENT  FORMULATION 


3.1  Geometry 


A  point  in  a  shell  structure  can  be  expressed  by  a  vector  sum  of  two  vectors.  One 
vector  is  a  position  vector  from  the  origin  of  the  coordinate  system  to  a  point  on  the 
reference  surface  of  the  shell  element,  and  the  other  vector  lies  from  the  end  of  the 
position  vector  to  the  point  under  consideration.  The  midsurface  of  the  shell  is  usually 
taken  as  the  reference  surface,  but  it  is  not  required  in  this  formulation.  The  second 
vector  described  above  is  usually  normal  to  the  reference  surface.  Therefore,  the  position 
vector  of  a  point  in  a  shell  can  be  expressed  as 

n  n 

=  (i=  1,2,3)  (1) 

k=l  ifc=l 

where  Xi  is  the  position  vector  of  a  generic  point  of  a  shell;  t],  and  C  are  the  natural 
coordinate  axes;  and  jy*(C)  are  two-dimensional  and  one-dimensional  shape  func¬ 

tions  in  the  natural  coordinate  system,  respectively;  xf  is  the  position  vector  of  the  node 
k  in  the  reference  surface;  is  the  unit  vector  at  the  node  fc;  and  n  is  the  number  of 
nodes  per  element.  In  the  present  formulation,  a  four-noded  shell  element  is  considered. 

Further,  the  unit  vector  is  defined  as 


_  ^^k^hottom 

j|  {xfy^p  —  (xf 


(2) 


where  top  and  bottom  indicate  the  top  and  bottom  surfaces  of  the  shell,  and  ||  ||  denotes 
the  Euclidean  norm.  The  one- dimensional  shape  function  H’‘  is  expressed  as 


hho  =  [J(i + 0(1  -  0  -  \ii  -  0(1 + 0]  II  II  (3) 


in  which  C  indicates  the  location  of  the  reference  surface  and  varied  from  -1  to  1.  C=0 
denotes  the  midsurface. 


3.2  Displacement 

The  displacement  field  in  a  shell  can  be  written  as 

.  n  n 

=  Y^N\trj)uf  +  +  VM,)  (i  =  1,2,3)  (4) 

h=l  *=1 
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Xi 

Ui 


Figure  1.  Shell  Element 

in  which  u,-  is  the  displacement  along  the  x,-  axis,  uf  is  the  nodal  displacement  at  the 
node  and  unit  vectors  and  14*-  lie  along  the  reference  surface.  and  V^i  are 

perpendicular  one  another.  and  0^^  are  rotational  degrees  of  freedom  along  the 

unit  vectors  Vi*-,  14*-,  and  14*-,  respectively.  The  right-hand  rule  is  assumed  for  the  positive 
direction  of  each  rotation.  That  is,  as  the  thumb  of  the  right-hand  is  in  the  direction  of 
each  unit  vector,  the  rotational  direction  of  the  hand  is  the  positive  rotation.  (See  Figure 
1-) 

and  0^^  are  the  bending  rotations  while  0^^  is  called  the  drilling  rotational  degree  of 
freedom.  To  explain  the  role  of  0^^  in  more  detail,  let  us  consider  a  flat  plate  whose  plane 
is  parallel  to  the  a;ia;2-plane.  Then,  the  displacement  field  in  Eq.  (4)  can  be  rewritten  as 

Ui  =  +  X30i  (i  -  1, 2, 3)  (5) 

where  and  M2  are  the  inplane  displacements,  and  M3  is  the  transverse  displacement. 
Superscript  mid  denotes  the  midplane  of  the  plate.  Equation  (5)  indicates  that  the 
transverse  displacement  is  not  constant  through  the  plate  thickness  (i.e.  along  the  x^- 
axis).  Therefore,  in  the  present  shell  formulation,  the  transverse  normal  strain  is  included 
as  well  as  the  transverse  shear  strains. 
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3.3  Strain-Displacement  Relation 

Six  strain  components  are  computed  from  Eq.  (4)  by  taking  derivative  with  respect 
to  the  ajj-axis.  The  result  can  be  written  in  matrix  form  hke 


{.}  =  [5]{d} 


(6) 


where 


{f}  —  {fll  f22  £33  7l2  723  7l3 


and 


{d}  =  {  rfi 

The  detailed  expression  for  [B^\  is 


5"] 


-  ajy* 

dxi 

0 

0 

-9^iV2\ 

9lV^i 

9lVl 

0 

dN’’ 

dx2 

0 

-9lV,% 

9lV,\ 

9lVi2 

0 

0 

dN’’ 

dxz 

-9lV^z 

9lV^ 

9lViz 

dN’’ 

dX2 

dN'’ 

dxi 

0 

1 

1 

9lV^i  +  9\Vh 

g^vi.  +  g^v^^ 

0 

dN” 

dSj 

dN’’ 

dX2 

-9lV^2-9lViz 

9lV,\  +  glV^^ 

glviz  +  g^v^ 

dN’’ 

-  dxs, 

0 

dN'’ 

dxi 

-9lv^i-gMz 

9lV^i  +  9lV^2 

glvi.^glv^ 

in  which 


In  addition,  the  vector  {d*}  is 


dxi 


dxi 


(7) 

(8) 

(9) 


(10) 


(11) 


{d*^}  —  {«!  U2  Uz  di  $2  dz} 


(12) 


3.4  Jacobian  Matrix 


In  order  to  compute  the  derivatives  such  as  and  it  requires  the  Jacobian 
matrix  defined  as 


where 


[j]  = 

Xz,( 

Xz,ti 

^2,c 

^3.C. 

dxi 

di  '' 

n 

=  E 

fc  =  l 

di 

n 

+E 

Jb=l 

dN^  ^ 
di  ^ 

dxi 

dr) 

n 

=  E 

k=zl 

dN^  i 
dr) 

n 

+E 

fc=l 

T» 

07] 

'Vi, 

(13) 

(*  =  1,2,3) 

(14) 

(*■  =  1,2,3) 

(15) 

4 


(16) 


Inverse  of  the  Jacobian  matrix  is  called  matrix  [i?].  Then, 

dN’^  „  „  dN’= 

dxi  di  (i- 1,2,3) 

dH’^  _  dH^ 

(e-1,2,3) 

Here  Rij  is  the  component  of  the  matrix  [R]. 


3.5  Element  Stiffness  Matrix 


Stresses  and  strains  can  be  expressed  as 

<Tij  =  -pSij  +  Sij  (i,  j  =  1, 2, 3)  (19) 

^ij  —  3^^*^  ^ij  (*iJ  —  1,2,3)  (20) 

where  p  and  Sij  are  the  hydrostatic  pressure  and  the  deviatoric  stresses,  and  e  and  e,j  are 
the  dilatation  and  the  deviatoric  strains. 

The  principle  of  virtual  work  states 

6W-8U  =  0  (21) 

where  6W  and  8U  are  the  virtual  work  of  the  external  and  internal  loads,  respectively.  In 
more  detail, 

8U  =  J  <Tij8eijdQ  (22) 


where  indicates  the  shell  volume. 


Substitution  of  the  stress  and  strain  expressions,  Eqs  (19)  and  (20),  into  the  virtual 
change  in  the  strain  energy,  Eq.  (22),  yields 


/  {-p8e  +  Sij8eij)d^ 
Jn 


For  linear  elastic  materials,  the  following  constitutive  equations  are  assumed. 


p  =  —Ke 


Sij  —  ^Gcij 
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where  K  is  the  bulk  modulus  and  G  is  the  shear  modulus.  Substitution  of  Eqs  (24)  and 
(25)  into  Eq.  (23)  results  in 


SU=  {Ke6e  +  2Geij6eij)dQ  (26) 

Jci 

Equation  (26)  is  derived  for  a  three-dimensional  solid.  For  a  shell,  plane  stress  con¬ 
dition  is  assumed  for  the  inplane  stress  components  o-n,  <722,  and  an-  With  this  imple¬ 
mentation,  Eq.  (26)  yields  the  element  stiffness  matrix  for  a  four-noded  shell  as  below: 

[  [Bpf[Dp][Bp]da+  l[Bf[D][B]dQ  (27) 

Jn 

where  the  first  and  second  integrals  are  related  to  the  volumetric  and  distortional  strain 
energies,  respectively.  Matrix  [B]  was  defined  in  Eq.  (8)  and  its  size  is  6  by  24  for  the 
four-noded  shell  element.  Matrix  [Bp]  is  the  upper  half  submatrix  of  the  matrix  [S],  which 
contains  the  matrix  components  related  to  three  normal  strains.  Thus,  matrix  [5p]  has 
the  size  of  3  by  24. 

In  addition,  matrices  [.D]  and  [Dp]  are  computed  from 


[D]  =  [Tf[D][T] 


(28) 


and 

[Dp]  =  [Tpf[Dp][Tp] 

in  which  the  transformation  matrix  [T]  is  expressed  as 


m  = 


(29) 


111 

ih 

^13 

/11/12 

/12/13 

/13/11 

^21 

^22 

^23 

/21/22 

h^hz 

hzhi 

^31 

^32 

^33 

/31/32 

hzhz 

/33/31 

2/11/21 

2/12/22 

^hzhz 

(/11/22  +  /21/12) 

(/12/23  +  /22/13) 

(/13/21  +  /23/11) 

2/21/31 

2/22/32 

^hzhz 

(/21/32  4-  /31/22) 

{hzhz  +  hzhz) 

(/23/31  +  /33/21) 

>2/31/11 

2/32/12 

‘^hzhz 

(/31/12  +  /11/32) 

{hzhz  +  /12/33) 

(/33/11  +  /13/31)  - 

(30) 


and  the  matrix  [Tp]  is  a  3  by  3  matrix  consisting  of  the  first  three  rows  and  columns  of 
the  matrix  [T].  Here,  lij  is  the  direction  cosines  of  the  unit  vector  Vi  with  respect  to  the 


Xj-axis. 


Matrix  [£>]  for  an  isotropic  elastic  material  is  defined  as 


[D]  =  2G 


p  (2Ga+2) 
3 

(2Go-l) 

3 

0 

0 

0 


(2Ga-l) 

3 

(2Ga+2) 

3 

0 

0 

0 


_  1 
3 

_i 

2^ 

3 

0 

0 

0 


0 
0 
0 

0 

0  0 


0 

0 

0 

0 

0 


(31) 
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where 


a 


+  ^ 
3  ^  2G 

2  I  ^ 

3  2G 


(32) 


and  K  is  the  transverse  shear  correction  factor.  The  other  matrial  property  matrix  [Dp] 
for  an  isotropic  elastic  material  is 


[^p]  = 


{\-a)K 
(1  -  a)K 


1 


{l-a)K  1 
(1  -  a)K  1 
1  1 


(33) 
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3.  EXPLICIT  TIME  INTEGRATION 


For  transient  analysis  of  a  shell  using  the  exphcit  time  integration  technique,  it  is  not 
necessary  to  form  element  stiffness  matrices  explicitly.  Instead,  element  internal  force 
vectors  are  computed  and  assembled.  Then,  the  acceleration  vector  is  computed  from 

{uy^[M]-\{Fe.ty-{Finty)  (34) 

where  {t/}  is  the  system  acceleration  vector,  [Af]  is  the  system  mass  matrix,  {Fe*t}  is  the 
system  external  force  vector,  is  the  system  internal  force  vector,  and  superscript  t 
denotes  time  at  t.  If  the  mass  matrix  is  a  diagonal  matrix,  the  invers  process  becomes 
very  simple  and  quick. 

The  system  internal  force  vector  is  obtained  from 

{Finty  =  (35) 

in  which  the  summation  is  over  the  entire  elements  and  {/}  is  an  internal  force  vector  for 
each  element.  The  force  vector  is  computed  from 


where  [Bp]  is  the  matrix  relating  deviatoric  strains  to  the  nodal  displacement/rotation 
vector,  and  {s}  is  the  deviatoric  stress  vector. 

When  the  constitutive  equation  for  voids,  such  as  the  Gurson  yield  function,  is  in¬ 
cluded  in  the  analysis,  the  pressure  in  the  above  equation  can  be  computed  from  Eq. 
(12)  in  Ref.  [1].  Because  Ref.  [1]  has  a  detailed  development  of  the  Gurson  model,  it  is 
not  repeated  here.  In  the  next  phase  of  research,  the  Gurson  model  will  be  implemented 
into  the  program. 


8 


4.  VERIFICATION  EXAMPLES 


Some  example  problems  were  solved  to  verify  the  present  shell  formulation.  Static 
analysis  was  conducted  for  a  plate  and  a  shell.  Static  analysis  checks  the  stiffness  matrix 
fomulation  which  is  the  essential  part  of  the  present  shell  element. 

The  first  example  was  a  clamped  square  plate  with  a  concentrated  force  at  the  center 
of  the  plate.  The  plate  dimension  is  10  in.  by  10  in.  with  the  thickness  of  0.1  in.  The 
elastic  modulus  and  Poisson’s  ratio  are  10  msi  and  0.2,  respectively.  The  applied  force  is 
40  lb.  The  analytical  deflection  at  the  center  of  the  plate  is  2.45x10-2  in.  For  the  finite 
element  modeling,  a  quarter  of  the  plate  was  considered  because  of  symmetry.  Using  four 
shell  elements  (2  elements  x  2  elements)  resulted  in  the  center  deflection  of  2.28x10-2  in. 
which  had  a  seven  percent  error  compared  to  the  analytical  solution.  On  the  other  hand, 
use  of  16  elements  (4  elements  x  4  elements)  gave  the  deflection  of  2.43x10-2  in.  which 
had  an  error  less  than  one  percent.  The  plate  was  modeled  to  be  in  the  xy-plane,  yz-plane 
or  xz-plane  to  check  the  transformation  matrix.  All  cases  gave  the  same  deflection  except 
that  its  direction  was  different  depending  on  the  orientation  of  the  plate  and  the  load. 


P 


R=5.0in.  t=0.094in.  L=10.35in. 
E=1.05e7psi  v=0.3125 


Figure  2.  Pinched  Shell  With  2x2  Mesh  Shown 
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The  next  problem  was  a  pinched  cyHnder  along  a  diameter  as  shwon  in  Figure  2. 
Both  ends  of  the  cyHnder  were  free.  The  cylinder  had  the  radius  of  5  in.,  the  length  of 
10.35  in.,  and  thickness  of  0.094  in.  The  elastic  modulus  and  Poisson’s  ratio  were  10.5  msi 
and  0.3125,  respectively.  The  pinched  load  was  100  lb.  Because  of  symmetry,  one  eighth 
of  the  cylinder  was  modeled.  The  inextensional  shell  theory  gave  the  radial  contraction 
of  0.1 11  Tin.  However,  the  solution  did  not  include  the  extensional  deformation.  As  a 
result,  the  finite  element  solution  including  extensional  deformation  was  expected  to  be 
higher  than  the  solution.  The  present  analysis  showed  that  the  four  element  model  (2 
elements  x  2  elements)  gave  the  displacement  0.0913  in.  As  the  mesh  was  refined,  the 
displacement  was  0.1219  in.  for  the  100  element  model  (10  elements  x  10  elements)  and 
0.1236  in.  for  the  256  element  model  (16  elements  x  16  elements),  respectively. 
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5.  CONCLUSIONS 


The  present  phase  of  this  research  was  to  formulate  a  shell  element  which  could  be 
used  to  model  the  void  growth  as  a  shell  structure  underwent  nonlinear  failure  process 
caused  by  a  shock  loading.  Gurson’s  model  was  chosen  to  simulate  the  void  growth  in 
a  metallic  material.  The  void  model  was  based  on  the  mean  stress  (or  pressure).  As  a 
result,  the  shell  element  formulation  was  developed  in  terms  of  separate  mean  stress  and 
deviatoric  stress  components  so  that  the  void  model  could  be  incorporated  into  the  shell 
element.  The  developed  shell  element  was  verified  using  some  analytical  problems. 

The  next  phase  of  the  work  is  to  further  verify  the  formulation  and  to  implement  the 
void  model  into  the  shell  element.  Once  this  phase  is  completed  successfully,  the  whole 
module  will  be  implemented  into  the  DYSMAS  program. 
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